Method and Apparatus for Evaluating Polynomials and Rational 



Functions 

Field of the Invention 

The present invention relates to a method and apparatus for computing values of 
5 mathematical expressions, and more specifically to a computer processing method and computer 
apparatus for computing binary representations of numerical values of polynomials and rational 
functions. 

.f ij Background of the Present Invention 

"Kara" 

Polynomials and rational functions, which are quotients of polynomials, are used, for 
yio example, by various branches of applied science for determining numerical values of 
N* mathematical expressions for modelling a property of a physical system, for example a rate of 

Jjjj air flow over an airfoil surface. Polynomials can be used to approximate complicated 

O 

y' mathematical expressions. Various methods, for example Horner's Rule that was disclosed in 
1819, are used for computing values of polynomials, and provide a degree of accuracy that may 
15 not be acceptable in certain situations, depending on the particular polynomial and the accuracy 
required. G.E. Forsythe et al in "Computer Methods for Mathematical Computations " 
Prentice-Hall (1977) discloses, in Section 4.2, using Horner's Rule for computing a numerical 
value of a polynomial. 

Steps for computing binary representations of numbers can create an unacceptably large 
20 deviation between a computed binary representation and its theoretical numerical value due to 
successive rounding errors. This can be an intolerable situation when a higher degree of 
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accuracy is required. Various methods can improve computation accuracy but they may require 
a significant increase in processing time and/or hardware. 

Agarwal et al in "New Scalar and Vector Elementary Functions for the IBM 
System/370 " IBM Journal of Research and Development: Vol.30 No.2 (March 1986), and Gal 
et al in " An Accurate Elementary Mathematical Library for the IEEE Floating Point Standard " 
ACM Transactions on Mathematical Software Vol. 17 No.l (March 1991) pp. 26 to 45, appear 
to disclose methods that include a lookup table having non-equidistantly spaced entries for 
computing values for elementary mathematical functions. 

Markstein in " Computation of Elementary Functions on the IBM RISC System/6000 
Processor " IBM Journal of Research and Development Vol.34 No.l (January 1990) appears to 
disclose a method for computing values of a function g(x) near a given point x h 

Gustavson et al in " Elementary Function Subroutines " Department of Mathematical 
Sciences of IBM Watson Research Center (January 1985) appears to disclose a program 
implementation for a software mathematical library. 

Object of the Present Invention 

An object of the present invention is to obviate the disadvantages of the prior art. 

Summary of the Present Invention 

A first aspect of the present invention provides a machine-processing method for 
computing a property of a mathematically modelled physical system, the steps including: 
reading, via a machine processing unit, input data including a value for each identified ordered 
coefficient of a first polynomial p(x) representing the property, the polynomial p(x) being 
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expressed as p(x) = • x J ) where./ = 0 to w, a value of a quantity x, a value of a 
predetermined x,-, and a value of a predetermined p(xj) correspondingly paired with said 
predetermined x t \ building, via the machine processing unit, a value of a second polynomial c(x) 
having ordered coefficients, the second polynomial c(x) being expressible as: 
cC*) = E(Cjt-x*) where k=0 to (n-l) so that the first polynomial p(x) is expressible as: 
p(x) =p(xi) + {x-X(} *c(x), including the steps of: determining, via the machine processing 
unit, a value for each ordered coefficient of the second polynomial c(x) by generating a plurality 
of machine processing unit signals to determine each of the ordered coefficients of the second 
polynomial c(x) from: C k = T,(P(k+\+j) • *{) where y = 0to(/j- 1 -A) ; determining, via the 
machine processing unit, a value of the second polynomial c(x) by generating a plurality of 
machine processing unit signals to determine: c(x) = X(C* •**) where £ = 0 to (« - 1) ; and 
constructing, via the machine processing unit, a value of the first polynomial p(x) by generating 
a plurality of machine processing unit signals to determined) =/*(*/) + {*-*/} • and 
outputting, via the machine-processing unit, the value of the first polynomial p(x) representing 
said property of the mathematically modelled physical system. 

Preferably the machine-implementable method is further adapted so that a difference 
between x and X/ is sufficiently small to achieve a desired accuracy of a final computation of the 
machine representation of a numerical value of the first polynomial p(x). 

Preferably the machine-implementable method is further adapted so that the step of 
reading the input data includes reading, via the machine processing unit, the input data from a 
machine-readable medium. 
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Preferably the machine-implementable method is further adapted so that the ordered 
coefficients of the second polynomial c(x) are computed from a mathematical expression 
derivable from: C k = Tt(P(k+x+j) ■ *J) where j = 0 to (n - 1 - A:). 

Preferably the machine-implementable method is further adapted so that the 
mathematical expression is a mathematical recurrence expression. 

Preferably the machine-implementable method is further adapted so that the 
mathematical recurrence expression is a forward mathematical recurrence expression. 

Preferably the machine-implementable method is further adapted so that the 
mathematical recurrence expression is a backward mathematical recurrence expression. 

Preferably the machine-implementable method is further adapted to implement the 
backward mathematical recurrence expression by including further steps for: equating, via the 
machine-processing unit, a value of a highest ordered coefficient of the second polynomial c(x) 
to a value of an identified highest ordered coefficient of the first polynomial p(x) by generating a 
plurality of machine processing unit signals to determine: = Pn\ and computing, via a 
machine processing unit, a value for each lower ordered coefficient of the second polynomial 
c(x) by generating a plurality of machine processing unit signals to determine: 
Cjfc-i = %i • Cjt + Pk where k=(n-l) to 1 . 

Preferably the machine-implementable method is further adapted so that the 
predetermined x y is selected from a set of predetermined values of x t . 

Preferably the machine-implementable method is further adapted so that the 
predetermined x t is a closest member of the set to the identified x. 
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Preferably the machine-implementable method is further adapted so that the step of 
determining a value of the second polynomial c(x) is computed by using Homer's Rule. 

Preferably the machine-implementable method is further adapted for determining a 
value of a denominator polynomial q(x) having identified ordered denominator coefficients, the 
denominator polynomial q(x) being expressible as: q(x) - 2(2; • x j ) where j = 0tom, including 
further steps of: adapting the input data to further include a value for each identified ordered 
denominator coefficient of the denominator polynomial q(x) 9 a value of a predetermined q(xj 
correspondingly paired with the predetermined x& and determining, via the machine processing 
unit, a value of another polynomial d(x) having ordered denominator coefficients, the another 
polynomial dfx) being expressible as: d(x) = 2 (At • x k ) where k = 0 to (m - 1 )so that the 
denominator polynomial q(x) is expressible as: q(x) = q{xi) + {x-x/} • d(x) 9 and a value for the 
denominator polynomial is resolved. 

Preferably the machine-implementable method is further adapted so that the first 
polynomial p(x) is a numerator polynomial p(x) 9 and p(x)-p(x) is computed, and p(xj) is not read. 

Preferably the machine-implementable method is further adapted for determining a 
value of a rational function r(x) being expressible as a quotient of the numerator polynomial p(x) 

p(x) 

and the denominator polynomial q(x) expressed as r(x)= including further steps of: 

adapting the input data to further including a value of a predetermined K x d correspondingly 
paired with the predetermined x,-; constructing, via the machine processing unit, a value of the 
rational function r(x) by generating a plurality of machine processing unit signals to determine: 
r\ ( wi i^^h*il^ p( x )-p( x i) 

^)=^i)-o— 555— ) h — 555 — . 
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Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to a Bessel function. 

Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to an error function (ERF). 

Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to a complementary error function (ERFC). 

Preferably the machine-implementable method is further adapted so that the rational 
function is an approximation to a log gamma function (LGAMMA). 

Preferably the machine-implementable method is further adapted so that all values are 
machine representations of some numerical value, the machine processing unit is a computer 
processing unit, and each machine representation is a binary representation operable with the 
computer processing unit, and the machine-readable medium is a computer-readable medium. 

A second aspect of the present invention provides a computer-readable program product 
having computer executable instructions for instructing a computer, the instructions tangibly 
embodying the machine-implementable method of the first aspect of the present invention. 

A third aspect of the present invention provides a computer-readable mathematical 
software routine library including computer executable instructions for instructing a computer, 
the instructions tangibly embodying the machine-implementable method of the first aspect of 
the present invention. 

Preferably, the computer-readable mathematical software routine library is further 
adapted so that the library is operatively associated with a software programming language. 
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A fourth aspect of the present invention provides a machine for computing a property of 
a mathematically modelled physical system, the steps including: reading, via a machine 
processing unit, input data including a value for each identified ordered coefficient of a first 
polynomial p(x) representing the property, the polynomial p(x) being expressed as 
p{x) = YSJPj *x j ) wherey = 0 to n 9 a value of a quantity x, a value of a predetermined x h and a 
value of a predetermined p(xi) correspondingly paired with the predetermined building, via 
the machine processing unit, a value of a second polynomial c(x) having ordered coefficients, 
the second polynomial c(x) being expressible as: c(x) = 2(Q * x*)where k = 0 to (n - 1 )so that 
the first polynomial p(x) is expressible as: p(x) =p(xd + {x-xt} * c(x) 9 including the steps of: 
determining, via the machine processing unit, a value for each ordered coefficient of the second 
polynomial c(x) by generating a plurality of machine processing unit signals to determine each 
of the ordered coefficients of the second polynomial c(x) from: 
Ck = H(P(k+\+j)*x f i ) where j = 0 to (n - 1 - A:); determining, via the machine processing unit, a 
value of the second polynomial c(x) by generating a plurality of machine processing unit signals 
to determine: c(x) = £(C* -x k ) where k= 0 to (« - 1); constructing, via the machine processing 
unit, a value of the first polynomial p(x) by generating a plurality of machine processing unit 
signals to determined) =/K*i) + {x-xi} -c(x); and outputting, via the machine-processing 
unit, the value of the first polynomial p(x) representing the property of the mathematically 
modelled physical system. 

Preferably, the machine is further adapted so that a difference between x and x, is 
sufficiently small to achieve a desired accuracy of a final computation of the machine 
representation of a numerical value of the first polynomial p(x). 
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Preferably, the machine is further adapted so that the means for reading the input data 
includes means for reading, via the machine processing unit, the input data from a 
machine-readable medium. 

Preferably, the machine is further adapted so that the ordered coefficients of the second 
polynomial c(x) are computed from a mathematical expression derivable from: 

Ck = XtfVitO •*!) where J" = 0 to (w - 1 - A) . 

Preferably, the machine is further adapted so that the mathematical expression is a 
mathematical recurrence expression. 

Preferably, the machine is further adapted so that the mathematical recurrence 
expression is a forward mathematical recurrence expression. 

Preferably, the machine is further adapted so that the mathematical recurrence 
expression is a backward mathematical recurrence expression. 

Preferably, the machine is further adapted to implement the backward mathematical 
recurrence expression by further including: means for equating, via the machine processing unit, 
a value of a highest ordered coefficient of the second polynomial c(x) to a value of an identified 
highest ordered coefficient of the first polynomial p(x) by generating a plurality of machine 
processing unit signals to determine: C n ~\ = P n \ and means for computing, via the machine 
processing unit, a value for each lower ordered coefficient of the second polynomial c(x) by 
generating a plurality of machine processing unit signals to determine: 
Cjfc-i = Xi -Ck + Pk where k = (n-\) to 1. 
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Preferably, the machine is further adapted so that the predetermined x, is selected from a 
set of predetermined values of x,. 

Preferably, the machine is further adapted so that the predetermined is a closest 
member of the set to the identified x. 

Preferably, the machine is further adapted for determining means for determining a 
value of the second polynomial c(x) is computed by using Horner's Rule. 

Preferably, the machine is further adapted for determining a value of a denominator 
polynomial q(x) having identified ordered denominator coefficients, the denominator 
polynomial qfx) being expressible as: q(x) = -*0 where / = 0 to m 9 including further steps 
of: adapting the input data to further include a value for each identified ordered denominator 
coefficient of the denominator polynomial qfx), and a value of a predetermined qfx) 
correspondingly paired with the predetermined x,; and determining, via the machine processing 
unit, a value of another polynomial d(x) having ordered denominator coefficients, the another 
polynomial dfx) being expressible as: d(x) = ^(D k *x k ) where lc= 0 to (m - l)so that the 
denominator polynomial qfx) is expressible as: q(x) = q(*i) + {*-*/}• d(x) 9 and a value for the 
denominator polynomial is resolved. 

Preferably, the machine is further adapted so that the first polynomial pfx) is a 
numerator polynomial pfx), and pfx)-pfx) is computed, and pfx) is not read. 

Preferably, the machine is further adapted for determining a value of a rational function 
rfx) being expressible as a quotient of the numerator polynomial pfx) and the denominator 

p(x) 

polynomial qfx) expressed as rfx)= "^j" > including further steps of: adapting the input data to 
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further including a value of a predetermined K x d correspondingly paired with the predetermined 
%i\ and constructing, via the machine processing unit, a value of the rational function r(x) by 
generating a plurality of machine processing unit signals to determine: 

fix) = rixd • (1 - g(x) ) + ^ . 

Preferably, the machine is further adapted so that the rational function is an 
approximation to a Bessel function. 

Preferably, the machine is further adapted so that the rational function is an 
approximation to an error function (ERF). 

Preferably, the machine is further adapted so that the rational function is an 
approximation to a complementary error function (ERFC). 

Preferably, the machine is further adapted so that the rational function is an 
approximation to a log gamma function (LGAMMA). 

Preferably, the machine is further adapted so that all values are machine representations 
of some numerical value, the machine processing unit is a computer processing unit, each 
machine representation is a binary representation operable with the computer processing unit, 
and the machine-readable medium is a computer-readable medium. 

A fifth aspect of the present invention provides a machine having a computer-readable 
program product having computer executable instructions for instructing a computer to embody 
the machine of the fourth aspect of the present invention. 
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A sixth aspect of the present invention provides a machine having a computer-readable 
mathematical software routine library including computer executable instructions for instructing 
a computer to embody the machine of the fourth aspect of the present invention. 

Preferably the computer-readable mathematical software routine library of the sixth 
aspect of the present invention is further adapted so that the library is operatively associated 
with a software programming language. 

Detailed Description of the Figures of the Present Invention 

To illustrate aspects of the present invention, the following figures are used, in which: 

Fig. 1 depicts a programming flow chart for computing a value of a polynomial 
p(x) in accordance with the present invention; 

Fig. 2 depicts a programming flow chart for computing a value of a rational 
function p(x)/q(x) in accordance with the present invention; and 

Fig. 3 depicts a programming flow chart for using a subroutine for computing a 
value of a Bessel function J 0 (x) in accordance with the present invention. 
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Detailed Description of a Preferred Embodiment of the Present Invention 

The present invention will be described with reference to an exemplary context of a 

computer processing method and apparatus for computing a binary representation of a numerical 

value of a polynomial p(x). 

5 For embodiments of the present invention, references to computing or calculating values 

or numbers will be understood to refer to programming instructions for instructing a computer 
to compute binary representations of numerical values of a mathematical expression or variable. 
It is understood that computing machinery can manipulate and transform binary representations 
||J of numerical values or numbers. 

M 

C3io A polynomial p(x) can be mathematically expressed as shown in the following equation: 

Sj P&) = Z Pj* j where/ = 0 to n in which p(x) = P 0 +P\ • x + P2 -x 2 + ... +P n • ityw. 

!! = i! 

» 1 



C9 Ordered coefficients of polynomial ^(x) are P 0 , i*„. Polynomial/?^ has n+1 ordered 

coefficients and is said to have degree n. A numerical value of polynomial p(x) can be computed 
15 by reading each coefficient of polynomial p(x) and an identified x and using steps in Equation 1. 

n(n+l) 

Following Equation 1 directly would require a computer processor to perform 2 
multiplications and n additions. To significantly reduce the processing time, Horner's Rule can 
be applied to Equation 1, in which polynomial p(x) can be expressed in Equation la, as shown 
below: 

20 p(x)=P 0 +x(Pi + x(P 2 +x(...(P„-i +x-P n )...))) Equ. 
la 
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By following Equation la, a program can perform n multiplications and n additions to 
compute a number for polynomial p(x) using less processing time than a program that follows 
Equation 1. 

A polynomial p(x) at x = JC/ can be expressed as follows: 

p(xd = 2 PjA where j = 0 to n Equ. 
2 

Therefore, a difference between Equation 1 and Equation 2 can be expressed as follows: 
p(x) - p(xd = 2 Pj(x J - xQ where j = 0 to « 

Expanding Equation 3 can provide the following expression: 

/>(*) = 2 - */) -2 [** • xp'" 1 "^] } where y = 0 to «, k = 0 to (/ - 1) Equ. 

4 

For Equation 4, in the first sum j = 0 to w ? and in the second sum £ = 0 to (j - 1 ). It is 
appreciated that Equation 4 can be expressed as follows: 

p{x) - p(xi) = (x-Xi) • 2 {(2 ^(Jt+i+y) ' ^/) • where £ = 0 to(n - l),y = 0 to(n -\-k) Equ. 
4a 

For Equation 4a, in the first sum k - 0 \o{n - 1) ? and in the second sum7 = 0 to(« - 1 - k) 
. By following Equation 4a, a value for polynomial p(x) can be computed by reading 
coefficients of the polynomial (P b P 2 , ... , P„) 9 and x, jc,-, and j^fo). Equation 4a can be expressed 
as a combination of Equation 5 and Equation 6, as follows: 
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p{x)-p(xi)= {x-Xi}-c(x) = {*-*/} -E (C k -x k ) where * = 0 to 
5 

C k = ZP{k^vA wherey = Oto(n-l-*) 
6 

Equation 5 includes an expression for a new polynomial cfxj having coefficients, from 
Q, Q„./; 3 that are not known a priori and therefore need to be determined. Equation 6 can be 
used for determining values of the coefficients of new polynomial c(x) by involving the 
coefficients of polynomial p(x) and an identified x h Optionally, other expressions for 
determining values for each coefficient of new polynomial c(x) can be derived from Equation 6. 

It is appreciated that Horner's Rule could be used to compute a value for polynomial 
c(x) after coefficients of polynomial c(x) have been computed by using Equation 6 or optionally 
from another expression that could be derived from Equation 6. By using Horner's Rule, 
Equation 5 can be expressed as follows: 

p(x) - p(xi) = {x-Xi} • {Co +x(Ci +x(C 2 +x(...(C„_ 2 +x ■ Cfl-i)...)))} Equ. 
5a 

Processing steps that follow Equation 5a begin within the innermost set of parentheses 
and proceed outwards to the outermost brackets. Calculated coefficients of new polynomial c(x) 
can subsequently be used in Equation 5a for computing a numerical value of polynomial p(x). 

Optionally, unknown coefficients of new polynomial c(x) can be computed by applying 
the principle of mathematical recurrence on Equation 6. For example, a forward mathematical 
recurrence can allow computation of coefficients of polynomial c(x) by beginning with the 
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lowest ordered coefficient, C 0 , and proceeding to the highest ordered coefficient, C (n . }) . 
Optionally, a backward mathematical recurrence can compute coefficients of polynomial c(x) by 
beginning with the highest ordered coefficient, C (n . }) and proceeding to the lowest ordered 
coefficient, C 0 . With each step of backwardly computing a coefficient of polynomial c(x), a next 
outer most bracketed term of Equation 5a can be advantageously computed. 

Optionally, a backward recurrence for determining each coefficient of new polynomial 
c(x) can be realized by using Equation 7 and Equation 8. Each subsequent lower-ordered 
unknown coefficient of new polynomial c(x) can be computed by using Equation 8. 

Cn-i = Pn Equ. 
7 

C k _i - Xi -Ck+ Pk where k=(n-l) to 1 Equ. 
8 

To compute a numerical value of a polynomial p(x) 9 a computer can read input data that 
includes a value for each identified ordered coefficient of the first polynomial p(x), a value of a 
quantity x, a value of a predetermined x h a value of a predetermined p(x) that is correspondingly 
paired with said predetermined x h An optional step can select a nearest x t that is closest to an 
identified x, in which the nearest jc, is selected from a set of predetermined x h and a nearest 
predetermined p(x) corresponds to the nearest x h Optionally, coefficients of polynomial c(x) can 
be computed by following Equation 7 and Equation 8. A value for polynomial p(x) can then be 
computed by following Equation 5a. The operations of computing the coefficients of 
polynomial c(x) in Equation 8 can be interleaved with computing the nested parenthesis in 
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Equation 5a so that each time a new C k is computed by Equation 8 the next outer nested 
parenthesis of Equation 5a is computed next. 

Provided p(x) - p(x$ does not have a multiple root at x = x h the relative round-off error 
in computing p(x) -p(x) with Equation 5 can be kept small by keeping x sufficiently close to x h 
This can be achieved by having available a number of sufficiently closely spaced x r -'s covering 
the range of x to be handled, and using the nearest x, to x. To obtain a reasonable result that 
would be satisfactory for a given computational purpose, a difference between x and x, can be 
sufficiently small enough to achieve a desired accuracy of a computed value of polynomial p(x). 

Optionally, accuracy can be further improved by choosing x s such that p(xj has extra 
accuracy beyond the precision of the floating-point number system of the computer. This can be 
achieved by selecting x, to be a floating-point number, by exhaustively searching from a starting 
desired value for x h until an x, is found for which pfxj is very close to a floating-point number 
(i.e. closer than half the distance between adjacent floating-point numbers). Such a 
floating-point approximation to p(x) has an accuracy beyond the precision of the floating-point 
number system. 

Fig. 1 shows the preferred embodiment of the present invention, which is represented as 
a computer processing method embodied in a programming flow chart for computing a machine 
representation of a numerical value of a polynomial p(x). Optionally, coefficients oip(x) can be 
placed in a vector of length n. The method optionally includes a step of accessing a table having 
a set of predetermined values of x r - and a corresponding set of predetermined values of p(x) each 
paired with a corresponding x,-. A table can contain values for the following expressions: 
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x h p{xt) for i=l to AT. 
8a 

In step 10, the method of the preferred embodiment begins. In step 15, each coefficient 
of polynomial p(x) and an identified x can be read by a computer processing unit. In step 20, the 
computer optionally selects a predetermined x t from a set of predetermined values of 
Preferably, the predetermined x t is a closest member of the set to the identified x. In step 30, the 
initial conditions for a programming loop are set. Variable ksum is an ongoing summation that 
represents the sum in Equation 5a. Variable ksum is initially set equal to coefficient P n of 
polynomial p(x). Variable ck is a temporary value for any coefficients of polynomial c(x) 9 in 
which variable ck is computed to be a specific coefficient of polynomial c(x) that is determined 
in step 50 as will be shown below. Variable ck is initially set equal to P m which is the highest 
ordered coefficient of polynomial p(x). Variable k is an index counter for a processing loop, in 
which variable k is initially set equal to (n-1), which is one less than the degree of polynomial 
p(x). In step 40, a query is performed to determine whether a value for variable k is equal to 
zero. If the numerical value for variable k is greater than zero, the method proceeds to step 50. If 
the numerical value for variable k is equal to zero, the method proceeds to step 60. In step 50, a 
value of variable ck is computed, which is coefficient C k by using the backward recurrence of 
Equation 8 that involves a higher ordered coefficient of polynomial p(x) 9 which is P h and a 
previously determined coefficient, C k . h of polynomial c(x). Next, an updated value for variable 
ksum is computed, which is an ongoing summation of the sum in Equation 5a, by determining a 
sum for the next outer bracketed term of Equation 5a. For remaining iterations of step 50, the 
summation represented in Equation 5a progresses from an innermost set of parentheses and 
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proceeds outwardly for each iterative step in an ongoing manner towards the outermost brackets. 
An updated value of variable k is computed so that the sum in Equation 5a can be resolved 
towards its outermost brackets. In step 60, the method uses a value for polynomial c(x) 9 which is 
expressed in Fig. 1 as variable ksum after variable k equals zero in step 40. A value of 
polynomial p(x) is computed as a rearrangement of Equation 5a, which is expressed as p(x) = 
p(Xj) + (x-Xi)*c(x). In step 70, the method for the preferred embodiment can stop or can be 
repeated from step 10. 

A second embodiment of the present invention provides a computer processing method 
for computing a machine representation of a numerical value of a rational function r(x). The 
rational function r(x) is expressed as a quotient of a numerator polynomial p(x) having a set of n 
predetermined coefficients from P 0 , P h P m so that p(x) = £ Pj-x! where/ = 0 to n, and a 
denominator polynomial q(x) having m predetermined coefficients from Q 0 , Q h Q m , so that 
q{x) = ^L Qj-x J where/ =0 to m . A rational function r(x) can be defined in the following 
expression: 

p(x) 

r ^=W) Equ ' 

8b 

From Equation 8b, it follows that: 

, . , W1 g(*H<*0 , P(x)-p(xj) 
r(x) = Kri) • (1 - ^ ) + 

The second method can optionally include access to a table having a set of 
predetermined values of x t and a corresponding set of predetermined values of q(x) each paired 
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with a corresponding x h and a corresponding set of predetermined values of r(x) each paired 
with a corresponding x h It is understood that each unique x g is paired with a correspondingly 
specific q(xi), and a correspondingly specific r(Xj). A spacing for table values of x, can be set 
such that a resulting absolute value of (x-Xj) is small enough to obtain a desired accuracy. The x t 
spacing required for a given desired accuracy is typically estimated with a perturbation analysis 
(using Taylor series, for example), and verified by numerical testing. The closeness of x to x,- 
needed to achieve a given accuracy in the computed r(x) depends on the particular rational 
function r(x). 

A table can contain values for the following expressions: 

xt 9 r(xi\ q(xi) for / = 1 to N Equ. 
10 

r{xi) and q(x{) for i-ltoN Equ. 
11 

Accuracy can be further improved by choosing the x's such that the following terms 
shown in Equation 11 can have extra accuracy beyond the precision of the floating-point 
number system of a computer. This can be achieved by selecting x, to be a floating-point 
number, by exhaustively searching from a starting desired value for x h until an x, is found for 
which so that both the quantities in Equation 11 are very close to floating-point numbers (closer 
than half the distance between adjacent floating-point numbers). Such floating-point 
approximations to r(xi) and q(x) have an accuracy beyond the precision of the floating-point 
number system. 
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The second method computes values of a numerator polynomial p(x) and a denominator 
polynomial q(x) by following Equation 5, as follows: 

^(x) -pipci) = {x -X;} • c(x) Equ. 
12a 

q(x) - q{xi) - {x-Xi} • d(x) Equ. 
12 

Polynomial c(x) has (n-1) unknown coefficients and polynomial dfx) has (m-1) 
unknown coefficients. The second method can determine coefficients of c(x) and d(x) and then 
determine a value for rational function r(x). 

Fig. 2 shows the second embodiment of the present invention, which is represented as a 
computer processing method embodied in a programming flow chart for computing a computer 
representation of a numerical value of a rational function r(x) = p(x)/q(x). Optionally, steps for 
accessing a table having a set of various predetermined values for x t a set of various 
predetermined values for q(x^ 9 and a set of various predetermined values for r(xj can be 
included. 

In step 110, the second method begins. In step 115, a computer reads values of each 
coefficient of polynomial p(x) and each coefficient of polynomial q(x) 9 and an identified x. In 
step 120, a predetermined x f from a set of predetermined values of x t can optionally be selected. 
Preferably, the predetermined x s is a closest member of the set to the identified x. These 
predetermined values can be optionally stored in a table or they can be supplied to the computer 
in some other means. Optionally, the table includes values for various xi, various q(xi) 
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associated with the corresponding x h and various K x d associated with the corresponding x t . 
Optionally, specific values for the x/s, q(xfs, and r(3^'s are chosen so that the qfxfs and r(xfs 
have extra accuracy beyond machine precision in the manner of the accurate table method that is 
disclosed by the prior art. It is appreciated that the computed result will be more accurate if the 
table has extra accuracy. Step 130, step 140, and step 150 are similar to step 30, step 40, and 
step 50 of Fig. 1, in that step 130, step 140, and step 150 provide a value for variable hum to be 
used in step 160 as will be shown next. In step 160, a value for variable dp (i.e., delta p) is 
computed,which is a computed value for p(x)-p(x^ for a numerator polynomial p(x). Variable dp 
will be used in step 205 as will be shown later. Step 170, step 180, and step 190 are similar to 
step 130, step 140, and step 150 of Fig. 2, in that step 170, step 180, and step 190 provide a 
value for variable ksum to be used in step 200 as will be shown next. In step 200 a value for 
variable dq (i.e., delta q), which is a computed value for q(x)-q(Xf) for a denominator polynomial 
q(x) is computed. Variable dq will be used in step 202 and 205 as shown later. In step 202 a 
numerical value for q(x) is computed by rearranging Equation 12. In step 205 a value for a 
rational function r(x) is computed by following Equation 9. In step 210 the second method can 
stop or can be restarted from step 110. 

A third embodiment of the present invention provides a computer processing method for 
computing binary representations of a value for a Bessel function. A programming flow chart 
can be developed by adapting Fig. 2 for computing a value for a Bessel function J 0 (x) 9 which is 
a Bessel function of the first kind of order 0. The third method can be further adapted to 
compute a value for various Bessel functions, for example J 0 , J h Y 0 , and Y h over a range of x 
(e.g. [0,8]) by first approximating the Bessel function by a piece-wise rational function. Bessel 
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functions can typically be approximated by rational functions having numerator and 
denominator polynomials with the same degree (i.e. same number of coefficients). 

To compute a value of a Bessel function J 0 (x) using the third method, the first task 
would be to determine a range that an identified x belongs in. The following is a list of steps for 
determining the proper range to which x belongs: if the identified x is a NAN (i.e., not a 

0 

legitimate number, such as q"), or +/- INF (i.e., any non-zero number divided by zero), the 
method does not compute a value for the Bessel function; if the identified x is less than zero, the 
method makes the sign of the identified x positive; if the identified x is greater than or equal to 
zero and less than 3.72* 10 9 , then the value of the Bessel function is 1, assuming the result is 
returned in IEEE double-precision floating point; if the identified x is greater than or equal to 
3.72* 10 9 and less than 1, then the method uses a Taylor series to compute a value for the Bessel 
function; if the identified x is greater than or equal to 1 and less than 1.5, then the method calls a 
subroutine that incorporates the present invention for determining a value for the Bessel 
function. The subroutine will be explained in Fig. 3. In each case where the subroutine is called, 
arguments P, Q, T, and n of the subroutine are read along with values of x and jc,- associated with 
the particular range that x belongs with; if the identified x is greater than or equal to 1 .5 and less 
than 1.9, then the method calls the subroutine for determining a value for the Bessel function; if 
the identified x is greater than or equal to 1.9 and less than 2, then the method calls the 
subroutine for determining a value for the Bessel function; if the identified x is greater than or 
equal to 2 and less than 4, then the method performs steps that are beyond the scope of the 
present invention. The steps involve accurate evaluation of a rational approximation to 
Jo(x)/(x-Zq), where z 0 is the first zero of J 0 \ if the identified x is greater than or equal to 4 and less 
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tlian 4.5, then the method calls the subroutine for determining a value for the Bessel function; if 
the identified x is greater than or equal to 4.5 and less than 5, then the method calls the 
subroutine for determining a value for the Bessel function; if the identified x is greater than or 
equal to 5 and less than 6, then the method performs steps that are beyond the scope of the 
present invention. The steps involve accurate evaluation of a rational approximation to 
J 0 (x)/(x-zj) , where z ; is the first zero of J 0 ; if the identified x is greater than or equal to 6 and 6.4, 
then the method calls the subroutine for determining a value for the Bessel function; if the 
identified x is grater than 6.4 and less than 7, then the method calls the subroutine for 
determining a value for the Bessel function; if the identified x is greater than or equal to 7 and 
less than 7.5, then the method calls the subroutine for determining a value for the Bessel 
function; if the identified x is greater than or equal to 7.5 and less than 8, then the method calls 
the subroutine for determining a value for the Bessel function; and if the identified x is greater 
than or equal to 8 then the method uses an asymptotic approximation for determining a value for 
the Bessel function, which is beyond the scope of the present invention. 

Fig. 3 shows the third embodiment of the present invention, which is represented as a 
computer processing method embodied in a programming flow chart for computing a machine 
representation of a numerical value of a Bessel function J 0 (x). 

In step 360 a subroutine, which can be called evalrat, of the flow chart of the third 
embodiment begins. The subroutine inputs an argument x, at which the value of a rational 
function is desired to be computed, and P and Q, which are vectors containing the coefficients 
of the numerator and denominator polynomials p(x) and q(x) 9 respectively, of the rational 
function, and 7, which is a table containing a set of predetermined x/s and corresponding r(x$ s 
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and qWs, and n 9 which is the degree of p(x) and q(x). In step 362 a computer processing unit 
reads values for each coefficient of numerator polynomial p(x) and denominator polynomial 
q(x) 9 in which both polynomials have the same number or coefficients, and also reads an 
identified x. In step 364, to achieve improved performance, a nearest x t and other required data 
5 are optionally read from a table. In step 366 various variables are initialized for subsequent use 
in an instruction loop. One software loop is used because the numerator and the denominator 
polynomials have the same degree (that is, share the same number of coefficients). This loop is 
similar to a combination of the two loops involving steps 130 to 150 and 170 to 190 in Fig. 2. In 
1^ step 368, a test determines whether the loop has been completed. In step 370, variables are 
©10 computed in a manner similar to that combining steps 150 and 190 in Fig. 1. In step 372, 

%£} 

^ variables are computed in a manner similar to that combining steps 160 and 200 to 205 shown 

-h 
\J 

l h i in Fig. 1. Step 374 marks the end of the subroutine of the flow chart of the third embodiment. 

For non-elementary special functions, for example the error function (ERF), the 
5 complementary error function (ERFC), or the Bessel functions,range-reduction properties, such 
|*15 as those available for the elementary functions, are not known. To directly apply prior art 
table-driven techniques, for example those used for the elementary functions, for these 
non-elementary special functions, would require a computer processing unit to read coefficients 
of each of a large number of polynomials, one associated with each table point. This can make 
the table unacceptably large. 

20 Computing steps that embody the methods of the present invention can be used to 

determine values for polynomials or rational functions approximating the non-elementary 
special function. Optionally, coefficients of a polynomial or rational function associated with 
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each individual table point, x h do not have to be stored; instead, these coefficients can be 
computed using simple formulas from only one (for polynomials) or two (for rational functions) 
stored function values per table point, resulting in tables of manageable size. In this sense, the 
present invention provides a benefit similar to that of range reduction formulas, for functions for 
which range reduction formulas do not exist. 

Mathematical software libraries, including those associated with a software 
programming language or operating system, can be adapted to instruct a general purpose digital 
computer to embody the present invention. It is appreciated that a software does not necessarily 
have to belong to a programming language or operating system. There are also stand-alone 
libraries of mathematical routines, for example, the IBM™ Engineering and Scientific 
Subroutine Library, ESSL™. 

There are many examples that illustrate the practical application of the present invention 
for mathematically modelling a property of a physical system. C. M. Close in ' The Analysis of 
Linear Circuits ", Harcourt, Brace & World (1966), discloses in Chapters 6 and 11 the use of 
polynomials for determining frequency response characteristics of electrical systems. The same 
concepts can be applied to non-electrical systems, such as mechanical, hydraulic, acoustical, and 
thermal systems. R.C. Weyrick in " Fundamentals of Automatic Control ", McGraw-Hill Book 
Company (1975), discloses in Chapters 2 and 6 the use of polynomials for mathematically 
modelling physical systems for predicting a behavior of a physical system. D. Halliday and R. 
Resnick in " Physics: Parts 1 and T \ John Wiley & Sons (1978) disclose in Chapters 3 and 4 the 
use of polynomials for mathematically modelling and predicting spatial co-ordinates of impact 
of a projectile being released from an aeroplane. 
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The present invention can be implemented in a computer-readable computer memory 
program having a medium for storing and carrying the program to a general purpose computer 
for instructing a computer processing unit to implement the embodiments of the present 
invention. A computer-memory product can be a floppy disk or compact disk that can be 
5 adapted to carry software that tangibly embodies the present invention. 



A computer processing method that tangibly embodies the present invention can be used 
for computing a number for a mathematical model and for comparing the computed number 
against an observed behaviour of a physical phenomenon. 



Q The concepts of the present invention can be further extended to a variety of other 

Oio applications that are clearly within the scope of this invention. 

%g Having thus described the present invention with respect to a preferred embodiment as 

;f implemented, it will be apparent to those skilled in the art that many modifications and 

enhancements are possible to the present invention without departing from the basic concepts as 
gjj described in the preferred embodiment of the present invention. Therefore, what is intended to 
15 be protected by way of letters patent should be limited only by the scope of the following 

claims. 
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